% X \in S^{n(p+1)}
% C \in S^{n(p+1)}
% TODO decompose into functions? not sure how to do that without losing cvx scope
% TODO vectorize; this works but it slow
function B = constrained_mlm(C, n, p, G)
	cvx_begin
		variable X(n*(p+1), n*(p+1)) symmetric

% construct D
		variable Y(n, n*(p+1));
		expression D(n, n*(p+1));
		for i=0:p
			D(1:n,1:n) = D(1:n,1:n) + X((1:n)+n*i,(1:n)+n*i);
		end
		for k=1:p
			for i=0:(p-k)
				D(1:n,(1:n)+n*k) = D(1:n,(1:n)+n*k) + 2*X((1:n)+n*i,(1:n)+n*(i+k));
			end
		end

% TODO why does this not work?
%		expression D(n, n*(p+1));
%		Y = compute_D(X, n, p);
%		D = Y;

% symmetrize/negate G and remove diagonal
		constrained_inds = repmat(logical(~(G+G') & (~eye(n))), 1, p); 

		minimize(-log_det(X(1:n,1:n)) + trace(C*X));
		X == semidefinite(n*(p+1));
		Y == D;
		Y(constrained_inds) == 0;
	cvx_end

	B0 = sqrtm(X(1:n,1:n));
	B1p = B0 \ X(1:n,(1+n):(n*(p+1)));
	B = [B0 B1p];
end
